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Establishing the phase diagram of hydrogen is a major challenge for experimental and the¬ 
oretical physics. Experiment alone cannot establish the atomic structure of solid hydrogen 
at high pressure, because hydrogen scatters X-rays only weakly. Instead our understanding 
of the atomic structure is largely based on density functional theory (DFT). By comparing 
Raman spectra for low-energy structures found in DFT searches with experimental spectra, 
candidate atomic structures have been identified for each experimentally observed phase. 
Unfortunately, DFT predicts a metallic structure to be energetically favoured at a broad 
range of pressures up to 400 GPa, where it is known experimentally that hydrogen is non- 
metallic. Here we show that more advanced theoretical methods (diffusion quantum Monte 
Carlo calculations) find the metallic structure to be uncompetitive, and predict a phase dia¬ 
gram in reasonable agreement with experiment. This greatly strengthens the claim that the 
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candidate atomic structures accurately model the experimentally observed phases. 

Hydrogen (H) is the simplest and most abundant of all elements and yet it displays amazing 
riehness in its phase behaviouit® it is observed to form a quantum erystalline state and orienta- 
tionally ordered moleeular phases, and it has been predieted to exhibit a liquid-metal phase at high 
pressures and low temperature^^, metallie superfluid and supereondueting superfluid state^^, 
and high-temperature superconductivit)^^^®. Several crystalline phases of solid molecular H have 
been observed in diamond anvil cell experiments carried out at pressures up to over 300 GPapStSI. 
The low-pressure phase I, which is a hexagonal close-packed structure formed of freely rotating 
molecules, transforms to a broken-symmetry phase II, in which the molecular rotations are re¬ 
stricted, at low temperature^^. The transition pressure decreases strongly with isotopic mas^^E^l^ 
and also depends on the total spin of the molecule^fEH. As the pressure is increased at low temper¬ 
atures, there is a further transition from phase II to a phase III at about 160 GPa, with the transition 
pressure for deuterium (D) exceeding that for H by about 12 GPcP^. Experimental studies have 
also demonstrated the existence of a phase IV at temperatures above a few hundred K and pres¬ 
sures above 220 GPali^GllliMiZI. Some constraints on the structures of the observed phases have 
been obtained from X-ray diffraction experiment^^^J^, but the low X-ray scattering cross section 
of H and the small sample sizes available limit the possible resolution. Infrared (IR) and particu¬ 
larly Raman spectroscopic measurements have yielded valuable information about the vibrational 
modes of H at high pressure^^^®^, but the available experimental data are insufficient to determine 
the structures of phases II, III, and IV. 
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Candidate structures for phases II, III, and IV have been suggested by structure searches 
based on density functional theory (DFXj^^^, although it should be emphasised that none of these 
structures has been identified as being unambiguously correct. The candidate structures for phase II 
consist of packings of molecules with bond lengths almost identical to the zero-pressure valuel^^lZl. 
We have modelled phase II using a molecular structure of Pljc symmetry with 24 atoms in the 
primitive unit cell, which we refer to as P2Jc-2A\ see Fig. [^a). (We adopt the convention of la¬ 
belling structures by their symmetry followed by the number of atoms per primitive cell.) P2\lc-2A 
is the most stable structure found to date in static-lattice DFT within the pressure range appropri¬ 
ate for phase II, and its vibrational characteristics are also compatible with those of phase II. We 
model phase III using a C2/c-24 structure consisting of layers of molecules whose bonds lie within 
the planes of the layers, forming a distorted hexagonal patteruP^, see Fig. [^b). This very stable 
structure can account for the high IR activity of phase IlP^. We also consider a molecular Cmca- 
12 structure!^, which is similar to C2/c-24, but slightly denser; see Fig. [^c). We model phase 
IV by a Pc-4S structureE^l^il, shown in Fig. [^d), which consists of alternate layers of strongly 
bonded molecules and weakly bonded graphene-like sheets. This type of structure was predicted 
by Pickard and Needd^. Pc-4S can account for the occurrence of stiff and soft vibronic modes in 
phase IV, and its stabilisation by temperature. Finally, we consider the Cmca-4 structure!^, which 
has weaker molecular bonds than C2/c-24 and Cmca-12, and is shown in Fig. [^e)- The main 
goals of our present work are to obtain accurate theoretical results for the relative stabilities of the 
P2i/c-24, C2lc-24, Cmca-12, Pc-4S, and Cmca-4 structures of H at pressures of 100-400 GPa and 
temperatures of 0-500 K, and to use these data to construct a temperature-pressure phase diagram 
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of H. We have not eonsidered phase I in our ealeulations, whieh is stable at low pressures, beeause 
an aeeurate deseription of this phase would require a full quantum treatment of the proton spin. 
Instead we foeus our attention on the phase behaviour at higher pressures, where the eandidate 
struetures are sueh that the nuelei are highly loealised and henee the motion of the protons is likely 
to be well-deseribed by eolleetive bosonie vibrational modes. 

Useful theoretieal deseriptions of solid H require very aeeurate ealeulations with an energy 
resolution of a few meV per atom. Various studies have shown that DFT eurrently eannot provide 
such accuracy for H structures, as evidenced by the disagreement of results obtained with differ¬ 
ent exchange-correlation functionals and the fact that DFT predicts H to be metallic at pressures 
above ~ 300 GPa, in contradiction with experimeni^HHElIlllllll. We have instead used the diffu¬ 
sion quantum Monte Carlo (DMC) methocP^to calculate static-lattice energy-volume relations for 
the different H phases. DMC is generally regarded as the most accurate first-principles method 
available for carrying out such studieP^®^. Furthermore, the low mass of the H atom means that a 
full treatment of quantum nuclear vibrational motion, including anharmonic effect^^5132l^ is crucial 
for an accurate description of the energetics. We have therefore used a DFT-based vibrational self- 
consistent field approachl^ to calculate anharmonic vibrational energies. We find that the use of 
DMC (and to a lesser extent the treatment of phonon anharmonicity) renders the metallic Cmca-A 
structure that is favoured in DFT energetically uncompetitive, leaving us with a phase diagram in 
reasonable quantitative agreement with experiment. 
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Figure 1: Atomic structures of the five H phases considered in this work, (a) P2\lc-2A, (h) 
C2/C-24, (c) Cmca-\2, (d) /’c-48, and (e) Cmca-A. The blue dumbbells indicate short bonds be¬ 
tween atoms (<0.8 A). The white dumbbells indicate long bonds between atoms (<0.9 A). The 
red lines indicate close contacts between atoms (<1.2 A) in the layered structures. P2\lc-2A con¬ 
sists of molecules arranged on a distorted hexagonal-close-packed lattice. C2/c-24, Cmca-\2, and 
Cmca-A consist of layers of molecules whose bonds lie within the planes of the layers, forming 
distorted hexagonal patterns, and we show top-down views of single layers, /’c-48 consists of 
alternate layers of isolated strongly bonded molecules and weakly bonded graphene-like sheets, 
and we show a top-down view of four layers. The structures are shown at a common DFT-PBE 
pressure (250 GPa). 
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Results 


Relative enthalpies Figure shows the static-lattice enthalpies of the structures relative to CUc- 
24. In Figs.j^a) and[^b) we report DFT enthalpies calculated using the Perdew-Burke-Emzerhof 
(PBEjiSl and Becke-Lee-Yang-Parr (BLYP) density functionald^HH. The relative DFT enthalpies 
are converged to better than 0.1 meV per atom with respect to k-point sampling and plane wave 
cutoff energy. The difference between the DFT-PBE and DFT-BLYP relative enthalpies arises 
chiefly from the energetics and not from the slightly different structures obtained from geometry 
optimisation calculations performed at fixed external pressures using the two different functionals: 
see Supplementary Note 1 and the accompanying Supplementary Fig. 1. In Fig. [^c) we report 
DMC enthalpies, which were obtained by fitting polynomials to the extrapolated infinite-system- 
size DMC energies as a function of volume, and differentiating the polynomials to obtain pressures. 
The structures used for the DMC calculations were obtained from DFT-PBE geometry optimisation 
calculations. We truncate the DMC enthalpy curves at the highest and lowest pressures at which 
we have performed calculations. 


The use of the DMC method has significant consequences for the static-lattice relative en¬ 
thalpies of the candidate structures. Compared with both DET-PBE and DET-BLYP, Cmca-4 and 
Cmca-\2 are destabilised with respect to C2/c-24, whereas P2\lc-2A is stabilised with respect to it. 


but in each case the DET-BLYP results are closer to the DMC enthalpies, as also found in Ref. 36 


Eor Cmca-4 and P2i/c-24 the difference between the DMC and the DET-BLYP results is greater 
than the difference between the DET-BLYP and DET-PBE results, while for Cmca-12 these differ- 
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Figure 2: DFT and DMC static-lattice enthalpy-pressure relations for the different H struc¬ 
tures relative to C2/c-24. (a) DFT-PBE, (h) DFT-BLYP, and (c) DMC. The relative DFT en¬ 
thalpies are converged to better than 0.1 meV per atom. The widths of the DMC lines indicate the 
estimated uncertainties in the enthalpies due to single-particle finite-size errors, which are greater 
than the uncertainties due to random sampling in the Monte Carlo algorithm, as explained in Sup¬ 
plementary Note 2. 
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ences are of similar size. Although DFT-BLYP happens to be relatively aeeurate in the pressure 
range of interest, it is elear that DFT is unable to provide a eonsistent, quantitative deseription of 
the relative enthalpies of the phases of H. 

Vibrational results The harmonie zero-point (ZP) eontributions to the enthalpies of the H phases 
inerease sublinearly with pressure, as shown in Fig.|^a), while the anharmonie eorreetions tend to 
deerease with pressure; see Fig. [^b). The harmonie ZP enthalpies are roughly thirty times larger 
than the anharmonie eorreetions. However, the differenees between the harmonie ZP energies of 
the five phases eonsidered at fixed pressure are similar in magnitude to the differenees between the 
anharmonie eorreetions, both being about 10 meV per atom, as shown in Figs.j^e) and[^d). This 
demonstrates that the variations in the anharmonie vibrational eorreetions are as important as those 
of the harmonie eontributions to the enthalpies in determining the relative stabilities of phases in 
this system. 

Structural phase transitions Figure shows the two struetural phase transitions that we have 
determined in this work, and our theoretieal temperature-pressure phase diagram for solid moleeu- 
lar H is shown in Fig. At 0 K, we find a transition from P2i/c-24 to C2/c-24 at around 235 ± 10 
GPa. The eorresponding transition pressure for D is 13 GPa higher. (Note that the differenee 
between H and D is purely due to the DFT vibrational free energy and henee the differenee in tran¬ 
sition pressures between H and D is relatively preeise.) Our transition pressure is around 75 GPa 
greater than those observed experimentally for the transition between phases II and III, but the 13 
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Figure 3: DFT-PBE vibrational contributions to the enthalpies of the H structures, (a) Har¬ 
monic ZP contributions to enthalpies, (b) anharmonic ZP corrections to enthalpies, (c) harmonic 
ZP enthalpies relative to C2/c-24, and (d) anharmonic ZP corrections relative to C2/c-24. P2i/c-24 
is destabilised by both harmonic vibrations and anharmonic corrections, relative to C2/c-24. Cmca- 
12, Cmca-4, and Pc-48 are all stabilised by harmonic vibrations but destabilised by anharmonic 
corrections, relative to C2/c-24. 
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Figure 4: Relative Gibbs free energies of the different H structures, (a) 0 K, (b) 150 K, (c) 

300 K, and (d) 400 K. The Gibbs free energies were ealeulated using statie-lattiee DMC ealeu- 
lations together with DFT-PBE harmonic and anharmonic vibrational calculations. The transition 
from PlJc-lA to CUc-lA occurs at around 235 ± 10 GPa between 0 K and 150 K. Pc-48 is sta¬ 
bilised by temperature with respect to C2/c-24. The complete set of relative enthalpies is shown in 
Supplementary Fig. 2. 
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Figure 5: Theoretical temperature-pressure phase diagram for H. The solid black lines show 
the phase transitions calculated in this work, i.e., the set of points at which the relative Gibbs free 
energy of two phases is zero. The dotted lines show the set of points at which the relative Gibbs 
free energy is one error bar from zero, and hence indicate the uncertainty in the phase boundaries. 
At pressures in excess of 350-375 GPa the Gibbs free energies of the C2/c-24 and /’c-48 structures 
are within error bars of each other. The grey region indicates the temperature-pressure conditions 
under which phase I is found to exist in experiments. 
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GPa difference between the transition pressures for H and D agrees well with the experimentally 
measured valueP^. We note that the theoretical transition pressures between H and D would only 
differ by around 6 GPa without the inclusion of anharmonic effects. 

As shown in Fig. Q we also find a temperature-driven transition from C2/c-24 to Pc-48 at 
pressures above 250 GPa and temperatures above 300 K, in good agreement with the experimen¬ 
tally observed transition between phases III and IV. In Fig. we show the relative free energies of 
C2/C-24 and Pc-48 at 300 and 400 K. At the lower temperature, C2/c-24 is marginally more stable, 
but at 400 K, Pc-48 has clearly become the more stable structure. The variation in the transition 
temperature with pressure is smaller than the uncertainty in that quantity, and so we report the 
C2/c-24-Pc-48 transition temperature as 320 ± 20 K. 

Discussion 

Our theoretical H phase diagram is in reasonable quantitative agreement with experiment, indicat¬ 
ing that the PlJc-lA, CUc-lA, and Pc-48 structures provide satisfactory models for phases II, III, 
and IV. These model structures reproduce the experimental Raman and IR spectra quite well. How¬ 
ever, there is a significant disagreement of about 75 GPa between the experimental and theoretical 
phase 11-111 transition pressure at 0 K. There are several possible reasons for our substantially 
larger phase II-III transition pressure. Firstly, the actual structure of phase III may be more sta¬ 
ble than the C2/c-24 model structure. However, C2/c-24 is the most stable non-metallic structure 
found in DFT searches over a wide range of pressures, and it is compatible with the Raman and IR 
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spectra of the observed phase III. If a significantly more stable structure than C2/c-24 were to be 
found for phase III, the excellent description of the transition from phase III to IV with increasing 
temperature obtained with our calculated data would be spoilt. Another possible explanation for 
the discrepancy with experiment regarding the phase 11-111 transition pressure could be that we 
have neglected a significant contribution to the energy of P2\tc-2A (our model for phase II). In 
particular, our calculations do not account for nuclear exchange effects, which are known to have a 
significant effect on the phase I-II transition pressure.^ However, reliable estimates of the size of 
nuclear exchange effects in solid H at high pressure are not currently available. Furthermore, nu¬ 
clear exchange effects are expected to be much smaller in D than H, because deuterons are bosons, 
whereas protons are fermions, and each deuteron has twice the mass of a proton. This suggests 
that nuclear exchange effects cannot be entirely responsible for the discrepancy in the phase 11-111 
transition pressure in both H and D. Our analysis of different finite-size corrections in Supple¬ 
mentary Note 2 (see also the accompanying Supplementary Figs. 3 and 4) indicates that finite-size 
effects in our relative enthalpies are well-controlled, but it is always possible that finite-size effects 
may be larger than anticipated. Finally, the fixed-node approximation is an uncontrolled source of 
error in our DMC calculations and, although fixed-node errors should largely cancel when relative 
energies are calculated, it cannot be ruled out that fixed-node errors may be larger in one phase 
than another. 

The results we obtain by combining our DMC static-lattice energies and harmonic and anhar- 
monic vibrational energies resolve a discrepancy between DFT and experiment for the transition 
between phases III and IV. The /’c-48 structure was proposed as a candidate for phase IV in 
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Refs. 28 and 29 because its Raman spectrum agrees well with the experimental one and because 


its weakly bonded layers lead to soft vibrational modes that thermally stabilise it. However, DFT 
static-lattice calculations together with the harmonic approximation for nuclear motion (used in 


Refs. 28 and 291 predict the Cmca-A structure to be energetically favoured at all temperatures in 


the relevant pressure range. (Note that Cmca-A is stabilised significantly by harmonic ZP energy; 
at the static-lattice DFT level it is not competitive, as shown in Fig. [^) The metallic nature of 
Cmca-A contradicts experiment, in which insulating structures containing strong molecular bonds 
are found up to pressures in excess of 300 GPa. The phase diagram predicted by DFT is shown in 
Supplementary Fig. 5. The use of static-lattice DMC energies and anharmonic vibrational energies 
destabilises Cmca-A, and we find that it is thermodynamically unstable over the entire pressure 
and temperature range considered here. Our results establish that DFT does not provide even a 
qualitatively correct description of the phase behaviour of hydrogen. We also find that Cmca-\2 
is unstable at the pressures and temperatures studied in this work. We have found an important 
discrepancy between our calculated phase TT-TTT transition pressure and experiment, which is cur¬ 
rently unresolved, although we have described possible physical reasons for the disagreement. Our 
calculations demonstrate that anharmonic vibrational effects are crucial for determining the relative 


stabilities of the phases. 


Methods 

Quantum Monte Carlo calculations The DMC methocP^Hlis capable of delivering much higher 
accuracy than DFT, and the scaling of the computational cost with system size enables the simula- 
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tion of the hundreds of atoms required for accurate calculations. We have used the DMC method to 
calculate static-lattice energies using H structures relaxed within DFT-PBE at a given external pres¬ 
sure. In DMC, the ground-state component of a trial wave function is projected out by simulating 
the Schrodinger equation in imaginary time, subject to the constraint that the nodal surface of the 
wave function is fixed to be that of the trial wave functiorl^^l^. We used Slater-Jastrow wave func¬ 
tions as implemented in the CASINO codeP^. Full technical details of our calculations can be found 
in Supplementary Note 2. The single-particle orbitals were obtained from the CASTEP codel^using 
the PBE exchange-correlation functional. The nuclei were represented by bare Coulomb potentials 
and appropriate cusp corrections were applied to the orbitals. We used a flexible Jastrow factoi® 
whose parameters were optimised using variational Monte Carlo (VMC j^. VMC and DMC simu¬ 
lations were performed using 96 and 768 atoms, and the results were extrapolated to infinite system 
size. Using the resources of the Oak Ridge Eeadership Computing Facility, we achieved statistical 
error bars of less than 0.3 meV per atom in all our DMC calculations. 

Anharmonic vibrational calculations We have calculated harmonic vibrational free energies 
by using the finite-displacement method to construct the matrix of force constants and diagonalis- 
ing the corresponding dynamical matrices over a fine vibrational Brillouin-zone grid, as described 
in Supplementary Note 3 (with accompanying data presented in Supplementary Figs. 6 and 7). 
We determined anharmonic corrections to the harmonic free energies using a vibrational self- 
consistent field metho(£^E3^, sampling the low-energy part of the DFT-PBE Born-Oppenheimer 
(BO) energy surface along harmonic normal modes to large amplitudes. The resulting anharmonic 
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Schrodinger equation for the nuclear motion was solved by expanding the wave function in a basis 
of simple harmonic oscillator eigenstates. Thermal occupation of excited states allowed us to cal¬ 
culate free energies at arbitrary temperatures. The vibrational free energy differences between the 
structures were converged to better than 1 meV per atom. Our approach does not describe possible 
melting. 
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